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Introduction 

Within the last few years several investigators 1 * 2 have established 
that in certain Mach number and angle of attack regions the conservative 
solution for steady full potential flow over a two-dimensional airfoil is 
nonunique. This phenomenon occurs at high subsonic Mach numbers and is 
connected with an indeterminateness of the shock location. Typically three 
distinct equilibrium flows exist, each with a different lift. For example, a 
symmetric airfoil at zero angle of attack might have a zero lift symmetric 
equilibrium, an asymmetric equilibrium with large positive lift, or its mirror 
image with negative lift. If the airfoil is placed at a small positive angle 
of attack one finds either a large positive lift, a small negative lift, or a 
large negative lift. None of these alternatives seems physically reasonable. 
Since the computed flow fields are essentially independent of the numerical 
scheme used, there is little doubt that the nonuniqueness is, in fact, a 
property of the differential equation. 

The present study deals with the evolution of these anomalous flows in 
time. In this context, nonuniqueness of the steady state arises from a 
sensitivity of the final equilibrium to changes in the initial conditions, or. 





equivalently, from the onset of instability in an equilibrium state. Previous 
related studies are reported in Ref. 3. 

The results given in Ref. 1-2 are for the steady full potential 
equation. The present work is based on the unsteady small disturbance 
approximation. This does change the equilibrium state structure of the 
problem somewhat, but, as will be apparent from the results, does not alter 
the basic multiplicity of the steady solutions. The problem is addressed by 
pure numerical experimentation using existing time dependent small disturbance 
algorithms. The intent of the paper is to report several new dynamic 
anomalies of unsteady small disturbance transonic flows, to verify the 
harmonic response results given in Ref. 3, (which were obtained with a 
slightly different algorithm) and, finally, to speculate on the physical 
relevance of these results. Attempts to improve the underlying physics of the 
mathematical model will be reported at a later date. 

In the first section of the paper the airfoil and flow parameter regime 
studied and the numerical algorithm employed are described. The airfoil is 
symmetric so that at zero angle of attack one expects the steady flow to be 
symmetric. These symmetric flows are described in the second section, and 
evidence is given that they are, in fact, weakly unstable over a range of free 
stream Mach numbers. 

The nature of the instability is mapped out in the third section. It is 
shown that, when the symmetric equilibrium is unstable, the flow evolves, very 
slowly, to either of two asymmetric equilibria, depending on the sign of the 
initial disturbance. The upper and lower surface shock waves, initially at 
the same location, split apart, one drifting downstream, the other upstream. 



The effect of this instability on dynamic response is then examined by 
imposing a harmonic pitch oscillation on the airfoil. The resulting 
hysteresis loops help to further clarify various anomalies which have been 
observed in the past (for example, non-zero mean lift in response to 
oscillation about zero mean angle^). 

In the fourth section, the effect of the boundary layer on the flow is 
examined. The results of a simple numerical test indicate that a fully 
interactive boundary layer does not significantly alter the stability of the 
inviscid equilibrium states. 

Description of the Flow Regime and Algorithm 

The problem examined in the present study is an NACA 0012 airfoil at and 
near zero angle of attack at Mach numbers between 0.8 and 0.9. This case was 
chosen to correspond to Salas' earlier work on steady flow. 2 

All the results reported herein were generated with the conservative, 
time accurate, small disturbance program XTRAN2L.4 developed at NASA 
Langley. This problem is a modification of LTRAN2-NLR 5 which includes 

(a) The <j>tt term in the wave equation, using the Rizzetta-Chin 
algorithm® 

(b) Non-reflecting far field boundary conditions? 

(c) A "monotone" difference switch to eliminate expansion 
shocks. 8 

The program was run with the NLR coefficients in the wave equation, 
without the Krupp scaling in the vertical direction, on an 80x40 grid (the 
default XTRAN3S xz grid). 

Since the present results are based on the small disturbance 
approximation, they are subject to transonic similarity rules. Hence most of 



the conclusions can be applied directly to any member of the OOXX series of 
airfoils. These scaling transformations are discussed in the Appendix. 
Symmetric Flow Results 

A baseline series of calculations was performed with the airfoil held 
fixed at zero angle of attack in a uniform free stream with Mach number 
ranging from 0.81 to 0.86. In the first run the airfoil was "turned on" in a 
uniform stream at M = 0.81 and the flow was allowed to develop to a steady 
state. This simulation spanned some 60 chords of travel (about 600 time 
steps). In the second run, the Mach number was re-set to 0.82, and the 
calculation continued from the 0.81 steady state until a new steady state was 
achieved. This process was repeated in Mach number increments of 0.01 to M = 
0.86. In each successive simulation, the flow appeared to reach a symmetric 
equilibrium in from ten to twenty chords of travel. The resulting steady 
pressure distributions are shown in Fig. 1. 

In reality, however, these calculations would not all converge to a 
symmetric equilibrium if continued indefinitely. If the apparently converged 
solution at M = 0.85, for example, is restarted with a large time step 
(corresponding to two chords of travel per step) the lift diverges 
exponentially, as shown in Fig. 2. In this case the lift history is, 
approximately. 


where Cl q is of order 10 -12 (since the initial asymmetry is a result of 
truncation error), and the growth rate, s, is roughly 0.01 U/c. These numbers 
explain why the instability was not observed in the original calculation: the 



initial flow asymmetry is exceedingly small and its time constant is large 
compared to the time scale for equilibration to the new Mach number. In fact 
it would take over 2000 chords of travel for lift coefficients on the order of 
0.001 to develop given the present initial conditions. 

It is important to note that the same growth rate is observed regardless 
of the time step (a large step is necessary only to capture a significant 
growth in a reasonable computational time). If the instability were numerical 
in origin, the growth rate would depend on At. We conclude that the 
instability must be a property of the differential equation. 

In order to examine the long time evolution of these unstable flows it is 
clearly desirable to introduce a controllable disturbance, rather than relying 
on truncation errors. 

Pulse Responses 

The stability of the symmetric equilibria described in the previous 
section was tested as follows: starting from the symmetric flow the airfoil 

was pitched up to 1/4° angle of attack and back to zero with a slow Gaussian 
pulse. The half width of the pulse was 40 and the time step 2 chords 
travelled - so that the pulse duration is comparable to the instability time 
scale and the time increment is small enough to resolve the disturbance and 
long enough to give the post-pulse history in a reasonable number of steps. 

The resulting dynamic lift responses are shown in Fig. 3 for Mach numbers 
0.82-0.86. It is clear from these results that the symmetric state is stable 
at M = 0.82, 0.83 and 0.86 - the disturbance generated by the pulse dies out 
in time, and the flow returns to its original symmetric zero lift condition. 

At Mach numbers between 0.83 and 0.86, though, the symmetric flow is 
unstable - the positive pulse in angle of attack induces a transition to an 
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asymmetric equilibrium state with positive lift (a negative pitch pulse would, 
of course, cause a transition to a mirror image negative lift state). The 
pressure distribution in the asymmetric equilibrium at M = 0.84, a = 0 is 

I 

shown in Fig. 4. Clearly the large lift is caused by a large shift in the 
upper/lower surface shock positions from their common symmetric state 
locations. Aside from the shock displacement there is very little change in 
the flow. 

Equilibrium lift and shock position results over the Mach number range of 
instability are shown in Figs. 5 and 6, respectively (each circled point was 
obtained by a pulse response calculation; starred points were taken from the 
initial symmetric flow fields). 

These figures show that at zero angle of attack there are three 
equilibrium flow fields at any Mach number between 0.835 and 0.858. The 
symmetric state is unstable and the two asymmetric states are stable (with 
respect to infinitesimal disturbances). Note further that the shocks cannot 
stand symmetrical ly anywhere between 72% and 88% chord. These limits are of 
interest because they are universal within the 00XX family of airfoils (x$/c 
is a valid transonic similarity parameter). It is important to observe as 
well that the asymmetric equilibrium states are not always associated with 
shocks standing near the trailing edge (as in Fig. 4). Indeed, near the lower 
bifurcation Mach number the stable shock positions are very close to the 
unstable symmetric location and well upstream from the trailing edge. Hence 
purely numerical problems associated with trailing edge shocks, while perhaps 
an issue, are not an explanation of the phenomenon. » 



Harmonic Responses 

One of the more important applications of small disturbance theory is the 
evaluation of the fluctuating aerodynamic forces which result from airfoil 
oscillation. Here we shall examine what happens to these forces when the 
airfoil is oscillated harmonically in pitch about zero angle of attack with a 
free stream Mach number in the range of instability. 

The harmonic response characteristics depend on the amplitude of the 
oscillation, the reduced frequency, k = oic/U, and, sometimes, on the way the 
oscillation is started. Figure 7 illustrates this, for the case M = 0.85, by 
means of dynamic Cl - a loops at three reduced frequencies; 0, 0.01 and 
0.05. 

The quasi -steady (k = 0) curve was obtained by a sequence of unsteady 
calculations (with different a) starting from the asymmetric a = 0 state 
labeled S in the figure. As long as a is greater than -0.06° the 
computation converges to a point on the upper branch of the quasi -steady 
curve. The lower branch is simply the mirror image of the upper, by 
symmetry. The middle branch, shown in dashes, is unstable and so cannot be 
computed with a time accurate algorithm. The quasi -steady curve is 
qualitatively similar to that obtained by Salas 2 , using a conservative full 
potential code. (He was also able to fill in the middle branch by using an 
inverse steady solver). 

A harmonic oscillation in a at an infinitesimal reduced frequency would 
lead to one of two results: 

(a) if the amplitude were less than 0.06° the flow would track the 

upper equilibrium curve (centered at S) or the lower curve, depending 
on initial conditions. 



(b) if the amplitude were larger than 0.06° the flow would 
transition between the upper and lower branches during 
the cycle, resulting in a hysteresis loop with abrupt 
jumps at a = jK).06°. 

Recall that the instability time scale is on the order of 100 c/U. Hence 
the quasi-static hysteresis loop just described would be seen only at reduced 
frequencies k « 0.01. The loop at k = 0.01 for a 1/4° amplitude 
oscillation is shown in Fig. 7. In this case the flow does transition back 
and forth between the two stable branches, though with considerable lag 
because the oscillation period is comparable with the instability time scale. 
Note that the mean lift is zero. At k = 0.05 (again for a 1/4° amplitude) 
the period is small compared to the instability scale. In consequence the 
flow does not have time to transition between branches during the cycle. The 
lift loop shown in Fig. 7 remains near the upper branch because the 
computation was begun on the upper branch. If the initial state had been 
chosen on the lower branch the final steady state oscillation would have been 
about the lower branch. Presumably any initial condition with k = 0.05 will 
eventually produce either the loop shown in Fig. 7 or its image. Note that at 
this "high" reduced frequency the mean lift is not zero. This type of result 
-in which a harmonic oscillation of a symmetric airfoil about zero mean angle 
produces a nonzero mean lift was observed and discussed by Dowell, et al. 3 
Boundary Layer Effect 

Experimental data 9 for the 0012 airfoil show no unequivocal evidence 
for any instability of the symmetric flow at a = 0 in the Mach number range 
where small disturbance theory predicts instability. One possible explanation 
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for this failure of the theory is that the interaction between the inviscid 
flow and the boundary layer stabilizes the shocks in the symmetric 
configuration. 

This hypothesis was tested with a preliminary viscous version of XTRAN2L 
currently being developed at NASA Langley (by James Howlett of the Unsteady 
Aerodynamics Branch). The model, which is based on Rizzetta's viscous 
modification of LTRAN2 10 , couples a quasi-steady, integral, turbulent 
boundary layer calculation of the displacement thickness with the existing 
inviscid code through the surface boundary condition. Strong interaction 
regions at the shock and trailing edge are not included in the formulation. 
From the point of view of the external flow, we see an inviscid flow over an 
airfoil with a compliant surface, such that the surface shape adjusts to the 
pressure distribution (and vice-versa). 

The test case chosen was this: 

1) The flow field was initialized in the inviscid asymmetric 
steady state at a = 0+, M = 0.85. 

2) The boundary layer was turned on and the flow field marched 
forward in time, simultaneously integrating the boundary 
layer and inviscid equations, with a and M fixed. 

Several outcomes of this test are conceivable: 

a) The flow approaches a symmetric equilibrium 

b) The flow oscillates without approaching an equilibrium 

c) The flow approaches a new asynmetric equilibrium. 

The actual outcome was (c): the boundary layer smears the shock profiles 

and displaces the shocks upstream, but does not alter the essential asymmetry 



of the flow. The initial and final pressure difference distributions are 
shown in Fig. 8. This shift in the loading corresponds to a drop in C|_ from 
0.307 with no boundary layer to 0.228 with a fully adjusted boundary layer. 

The transition to the new viscous equilibrium was essentially complete 
within ten chords of travel. As a precaution the calculation was continued to 
over a hundred chords with no significant change. 

Although the particular boundary layer model used in this calculation has 
several theoretical shortcomings, it is difficult to imagine that a more 
refined model would alter the basic conclusion: including boundary layer 

interactions in the small disturbance inviscid model does not cure its 
"unacceptable" dynamic properties. 

Conclusions 

The major conclusion to be drawn from the results presented here is 
clear: when solved as a weak conservation law, the transonic small 

disturbance equation predicts flow dynamics which are physically unexpected. 

At lower supercritical Mach numbers there is one, stable, equilibrium state. 

As the Mach number increases the original state becomes unstable and two new 
stable states appear. At a certain high Mach number these states 
re-coalesce. This sequence appears to be unaffected by (at least a simple 
version of) inviscid flow/boundary layer interaction. Moreover the 
equilibrium state structure (and presumably the dynamic response 
characteristics as well) is common to small disturbance and full potential 
theory. 

Perhaps the most surprising property of the instability is its extreme 
slowness - even on the time scale of upstream wave propagation. No 



satisfactory explanation for this long time behavior has been found, though it 
is clearly associated with a slow drift in the shock position. 

The most pressing unresolved question is whether these anomalous flows 
are intrinsic to inviscid theory or simply a consequence of some approximation 
made in small disturbance theory (and, presumably, in full potential theory as 
well). In this regard it is significant that both the small disturbance and 
full potential theories display the anomalous behavior while Salas' solutions 
of the Euler equations show no conclusive evidence of it. Small disturbance 
and full potential theory share a common neglect of the vorticity and entropy 
which are generated by shock waves - effects which are present in the Euler 
equations. Hence there is a suggestion that the shock instabilities are 
caused by the neglect of these mechanisms. 

A more persuasive argument can be made, however, that the instabilities 
must arise in the Euler equations as well. Small disturbance theory (unlike 
full potential theory which is ad-hoc when shocks are present) is a rational 
asymptotic limit of the Euler equations for Mach numbers near one and small 
thickness ratios (6 + 0, (1 - )/ ^2/3 fixed). In the transonic scaling 

the Euler equations can be written symbolically as 

N q (U;K) = 6^ /3 N 1 (U) 

where U is the scaled state vector, K is the transonic similarity parameter 
and N 0 and Ni are nonlinear operators. If we set 6 = 0 we get small 
disturbance theory (N 0 = 0). It is difficult to imagine how the solution of 
the Euler equations for 6 nonzero but arbitrarily small could differ 
substantially from the solution of small disturbance theory. 

If the solution of the Euler equations does approach the solution of 
small disturbance theory as 5 + 0, then the Euler equations must display the 



same instability and nonuniqueness properties as seen in small disturbance 
theory. Since Salas' numerical results indicate that the Euler solutions are 
well behaved for the 0006 and 0012 airfoils it would seem that a bifurcation 
must occur at some finite 6 less than 0.06 (for the NACA 00XX family at zero 
angle of attack). This might occur as sketched in Fig. 9. Figure 9a shows a 
hypothetical steady state bifurcation for a fixed K (between the limits 
imposed by small disturbance theory). Figure 9b shows a hypothetical M-S 
parameter map of the region of nonuniqueness in the steady Euler equations. 
These figures are consistent both with the present numerical results and with 
Salas' Euler solutions. Whether or not they in fact represent the behavior of 
the Euler equations, of course, awaits further investigation. If it does turn 
out that Figure 9 is qualitatively correct, then either the phenomenon is real 
or it results from the neglect of viscous forces. In the latter event, a 
better model of viscous effects would be required than that used in the 
present investigation. In this regard it should be noted that existing 
experimental evidence (which is for relatively thick airfoils) does not 
preclude the phenomenon illustrated in Fig. 9. It is quite possible that the 
flow instabilities reported here do in fact occur in nature, but that their 
domain of occurrence has simply not been studied experimentally. 

Appendix 

Similarity Laws for the Instability Region 
Since the instability is extraordinarily weak, it should be describable, 
with negligible error, by the low frequency small disturbance equation, 

2M<j> xt + (M 2 - 1 + B<t> x )<t» xx - 4> yy « 0 (A.l) 

where x,y are measured in units of the chord c, t in units of c/U, and where B 
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is a constant. Correspondingly we impose the quasi-steady boundary conditions 

<i> = xf'(x) on airfoil, 0<x<l, y = 0 (A. 2) 

a<|> =0 on wake , x > 1, y = 0 (A. 3) 

where t is a scale factor for both thickness and angle of attack. 

If no external time scale is introduced through the boundary condition 
(eg. if the airfoil is fixed) then the solution of the above problem must 
scale as: 


* "j? (x, y. t; K) 

3 = / 1 - f? , K = Bt/3 3 (A. 4) 

y = 3y, t * 3 2 t/M, 

from which it is apparent that the lift and shock position are: 


C L = J C L (t; K) 
x$ = x^(t"; K) 


(A. 5) 


Suppose, now, that a symmetric airfoil (thickness/chord ratio t) is at 
zero angle of attack. Suppose, further, that the flow is initially in the 
symmetric equilibrium state, with Cl = 0 and x$ = *s 0 (K). ^ this 

flow is unstable, so that the lift grows like e st , then the growth rate of 
the instability must scale like: 


S = K s(K) . (A. 6) 

3 

Instability onset occurs when ? = 0, that is, for specific values of the 
similarity parameter K, or, equivalently, at specific values of the symmetric 
shock position x$ Q . For the 0012 airfoil, with the NLR coefficient 
B = M 2 (2 + (y - 1)M 2 ) instability appears at M = 0.835 and disappears at 
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0.858. These limits correspond to Kj = 1.144, x$ 0 = 0.72 and 
l <2 = 1.496, x$ Q = 0.88 respectively. Hence, for any airfoil 00x, with 
any coefficient B, the symmetric flow will be unstable whenever Bx/g3 is 
between Kj and K 2 , or, equivalently, whenever the symmetric shock lies 
between 0.72 and 0.88. When the instability occurs the flow will transition 
to a new asymmetric equilibrium, the lift and shock positions of which can 
easily be calculated from (A. 5) and Figures 5 and 6. 

A puzzling feature of the numerical results is the extremely slow growth 
rate of the instability, which has a time scale t ~ 100. This is not 
explained by the scaling law since in the middle of the unstable region 
(M = 0.85), the scaled time is T ~ 33. We give here a simple ad-hoc model 
of the instability, based on the scaling law and the computed equilibrium 
shock locus (Fig. 6), which does predict the correct growth rates. 

Suppose that the shock velocity is governed by the differential equation 

dx 2 

lit = " A ~M ( X S ' X S “ b >< x S " X S ^ X S " X S + b ) ( A * 7 ) 

0 0 0 

so that x$ Q is the symmetric equilibrium position and x$ 0 ± b are the 
asymmetric equilibria. By similitude, b and A must be determined by x$ Q . 

Note that if A is positive the symmetric state is unstable and the two 
asymmetric states are stable. Linearizing about the symmetric state we find 
the behavior, 

x- - x- + Ce st (A. 8) 

a 0 

where 

s = AbV/M (A. 9) 
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From Fig. 6 the maximum value of b is 0.15 at M = 0.845, which yields a 
maximum growth rate of 0.008A. This value agrees with the observed growth 
rate if A is of order 1. Over the range of instability, then, we expect 
growth rates given by the empirical relation 

e 2 

S = 2.5 (Xc - .72) (.88 - Xo ) (A.10) 

n ^0 ^0 

(which comes from fitting a parabola for b 2 to the results of Fig. 6). 

This derivation is, of course, quite crude since Eq. (A. 7) has no 
justification beyond its consistency with Fig. 6 and the scaling laws. What 
has been shown is that the small growth rates observed arise "naturally," 
without large nondimensional factors being inserted arbitrarily. A rational 
analysis of the growth rate would require solving the eigenvalue problem which 
results from linearizing Eqs. (A. 1 )- (A. 3) about the symmetric equilibrium. 

Even if this were done, though, the cause of the instability would likely 
remain as mysterious as it is now. 
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(a) Lift bifurcation at fixed K = (1 -M^)/a. 


Figure 9. 



(b) Region of nonuniqueness. 
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